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Abstract: A scheme for rapidly and accurately computing solutions to boundary 
integral equations (BIEs) on rotationally symmetric surfaces in K 3 is presented. The 
scheme uses the Fourier transform to reduce the original BIE defined on a surface 
to a sequence of BIEs defined on a generating curve for the surface. It can handle 
loads that are not necessarily rotationally symmetric. Nystrom discretization is used 
to discretize the BIEs on the generating curve. The quadrature used is a high-order 
Gaussian rule that is modified near the diagonal to retain high-order accuracy for 
singular kernels. The reduction in dimensionality, along with the use of high-order 
accurate quadratures, leads to small linear systems that can be inverted directly via, 
e.g., Gaussian elimination. This makes the scheme particularly fast in environments 
involving multiple right hand sides. It is demonstrated that for BIEs associated 
with Laplace's equation, the kernel in the reduced equations can be evaluated very 
rapidly by exploiting recursion relations for Legendre functions. Numerical examples 
illustrate the performance of the scheme; in particular, it is demonstrated that for a 
BIE associated with Laplace's equation on a surface discretized using 320 000 points, 
the set-up phase of the algorithm takes 2 minutes on a standard desktop, and then 
solves can be executed in 0.5 seconds. 



1. Introduction 



This paper presents a numerical technique for solving boundary integral equations (BIEs) 
defined on axisymmetric surfaces in M 3 . Specifically, we consider second kind Fredholm 
equations of the form 



(1.1) 



a(x) + J k(x, x') a{x') dA(x') = f(x) 



xeF, 



under two assumptions: First, that T is a surface in M 3 obtained by rotating a curve 7 
about an axis. Second, that the kernel k is invariant under rotation about the symmetry 
axis in the sense that 

(1.2) k(x, x') = k{0 — 9', r, z, r', z'), 

where (r, z, 0) and (r', z' , 9') are cylindrical coordinates for x and x' , respectively, 

x = (r cos 9, r sin 9, z), 

x = (r cos 9', r sin9 , z), 

see Figure [lT] Under these assumptions, the equation ( |1.1[ ), which is defined on the two- 
dimensional surface T, can via a Fourier transform in the azimuthal variable be recast as 
a sequence of equations defined on the one-dimensional curve 7. To be precise, letting a n , 



f n , and k n denote the Fourier coefficients of a, /, and k, respectively (so that (2.6), (2.7) 
and (2.8) hold), the equation (1.1) is equivalent to the sequence of equations 



(1.3) o- n (r,z) + V2iT k n (r,z,r',z')o- n (r',z')r'dl(r' 1 z') = f n (r,z), (r, z) £ 7, n <E Z. 
Whenever / can be represented with a moderate number of Fourier modes, the formula 



(1.3) provides an efficient technique for computing the corresponding modes of a. The 
conversion of (1.1) to (1.3) appears in, e.g., [TS], and is described in detail in Section[2~l 
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Equations of the type ( 1.1 ) arise in many areas of mathematical physics and engineering, 
commonly as reformulations of elliptic partial differential equations. Advantages of a BIE 
approach include a reduction in dimensionality, often a radical improvement in the con- 
ditioning of the mathematical equation to be solved, a natural way of handling problems 
defined on exterior domains, and a relative ease in implementing high-order discretization 
schemes, see, e.g., [2]. 

The numerical solution of BIEs such as ( 1.1 ) poses certain difficulties, the foremost being 
that the discretizations generally involve dense matrices. Until the 1980s, this issue often 
times made it prohibitively expensive to use BIE formulations as numerical tools. However, 
with the advent of "fast" algorithms (the Fast Multipole Method O [10] , panel clustering 
|13j . etc.) for matrix- vector multiplication and the inversion of dense matrices arising 
from the discretization of BIE operators, these problems have largely been overcome for 
problems in two dimensions. This is not necessarily the case in three dimensions; issues such 
as surface representation and the construction of quadrature rules in a three dimensional 
environment still pose unresolved questions. The point of recasting the single BIE (1.1) 
defined on a surface as the sequence of BIEs (1.3) defined on a curve is in part to avoid 
these difficulties in discretizing surfaces, and in part to exploit the exceptionally high speed 
of the Fast Fourier Transform (FFT). 

The reduction of (1.1) to (1.3) is only applicable when the geometry of the boundary is 
axisymmetric, but presents no such restriction in regard to the boundary load. Formulations 
of this kind have been known for a long time, and have been applied to problems in stress 
analysis [3], scattering [HI EH I22J [23] , and potential theory [T21 EH QH [20] . Most of these 
approaches have relied on collocation or Galerkin discretizations and have generally relied 
on low-order accurate discretizations. A complication of the axisymmetric formulation is 
the need to determine the kernels k n for a large number of Fourier modes n, since direct 
integration of (1.2) through the azimuthal variable tends to be prohibitively expensive. 
When k is smooth, this calculation can rapidly be accomplished using the FFT, but when 
k is near-singular, other techniques are required (quadrature, local refinement, etc.) that 
can lead to significant slowdown in the construction of the linear systems. 

The technique described in this paper improves upon previous work in terms of both 
accuracy and speed. The gain in accuracy is attained by constructing a high-order quadra- 
ture scheme for kernels with integrable singularities. This quadrature is obtained by locally 
modifying a Guassian quadrature scheme, in a manner similar to that of [HE]. Numerical 
experiments indicate that for simple surfaces, a relative accuracy of 10 -10 is obtained using 
as few as a hundred points along the generating curve. The rapid convergence of the dis- 
cretization leads to linear systems of small size that can be solved directly via, e.g., Gaussian 
elimination, making the algorithm particularly effective in environments involving multiple 
right hand sides and when the linear system is ill-conditioned. To describe the asymptotic 
complexity of the method, we need to introduce some notation. We let Np denote the 
number of panels used to discretize the generating curve 7, we let Nq denote the number of 
Gaussian points in each panel, and we let Np denote the number of Fourier modes included 
in the calculation. Splitting the computational cost into a "set-up" cost that needs to be 
incurred only once for a given geometry and given discretization parameters, and a "solve" 
cost representing the time required to process each right hand side, we have 



(1.4) T setU p~ N$N%N F log(N F ) + N P N^N$ + Np N G Np 

construction of linear systems inversion of systems 
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and 

(1.5) T solve ~N P N G N F Iog(JV F ) + N$N G N F • 

V v ' V v ' 

FFT of boundary data application of inverses 



The technique described gets particularly efficient for problems of the form (1.1) in 
which the kernel k is either the single or the double layer potential associated with Laplace's 
equation. We demonstrate that for such problems, it is possible to exploit recursion relations 



for Legendre functions to very rapidly construct the Fourier coefficients k n in (1.3). This 



reduces the computational complexity of the setup (which requires the construction of a 



sequence of dense matrices) from (1.4) to 

T setup ~ iVjp N G N F log(iV F ) + iV P N G iV F + N$, N G N F . 

Numerical experiments demonstrate that for a problem with N F = 80, N G = 10, and 
iV F = 400 (for a total of 80 x 10 x 400 = 320 000 degrees of freedom), this accelerated scheme 
requires only 2.2 minutes for the setup, and 0.46 seconds for each solve when implemented 
on a standard desktop PC. 

The technique described in this paper can be accelerated further by combining it with 



a fast solver applied to each of the equations in (1.3), such as those based on the Fast 
Multipole Method, or the fast direct solver of [TTj. This would result in a highly accurate 
scheme with near optimal complexity. 



The paper is organized as follows: Section [2] describes the reduction of (1.1) to (1.3) 



and quantifies the error incurred by truncating the Fourier series. Section [3 presents the 
Nystrom discretization of the reduced equations using high-order quadrature applicable to 
kernels with integrable singularities, and the construction of the resulting linear systems. 
Section [4] summarizes the algorithm for the numerical solution of (1.3) and describes its 



computational costs. Section [5] presents the application of the algorithm for BIE formula- 
tions of Laplace's equation and describes the rapid calculation of k n in this setting. Section 
[6] presents numerical examples applied to problems from potential theory, and Section [7] 
gives conclusions and possible extensions and generalizations. 

2. Fourier representation of BIE 

2.1. Problem formulation. Suppose that T is a surface in M 3 obtained by rotating a 
smooth contour 7 about a fixed axis and consider the boundary integral equation 

(2.1) a(x) + J k(x,x')a(x')dA(x') = /(as), x e r. 

In this section, we will demonstrate that if the kernel k is rotationally symmetric in a sense 



to be made precise, then by taking the Fourier transform in the azimuthal variable, (2.1 ) can 
be recast as a sequence of BIEs defined on the curve 7. To this end, we introduce a Cartesian 
coordinate system in R 3 with the third coordinate axis being the axis of symmetry. Then 
cylindrical coordinates (r, z, 9) are defined such that 

x\ = r cos 0, 
X2 = r sin 9, 
x 3 = z. 



Figure 1.1 illustrates the coordinate system. 



The kernel k in (2.1) is now rotationally symmetric if for any two points x,x' £ T, 
(2.2) k(x,x') = k{9-9',r,z,r',z'), 

where (9',r',z') are the cylindrical coordinates of x' . 



1 




2.2. Separation of variables. We define for n £ Z the functions f n , a n , and k n via 

— in8 



(2.3) 
(2.4) 
(2.5) 



U(r,z) 
cr n (r,z) 
k n (r,z,r',z') 



T v27T 
—ind 



f(e,r,z)de, 



T v27T 

D —in9 



a(6,r,z)d9, 



T \/27r 



&;(#, r, z, r , z') <i#. 



The definitions (2.3), (2.4), and (2.5) define /„, a n , and k n as the coefficients in the Fourier 
series of the functions /, a, and k about the azimuthal variable, 

(2.6) 
(2.7) 
(2.8) 



tTi, V 2?r 



fc(a;, a;') = - 0', r, z, r', z') = ^ 



jn(e-e') 



2n 



k n (r,z,r',z'). 



To determine the Fourier representation of (2.1 ), we multiply the equation by e in8 /y/2n 
and integrate over T (for our purposes, we can think of T as simply the interval [—it, it]). 



Equation (2.1) can then be said to be equivalent to the sequence of equations 

—ind 



(2.9) 



7XT 



T v 2?r 



k(x,x')d9 



a(x')dA(x') = f n (r,z) 



n £ 



Invoking ( |2.8[ ), we evaluate the bracketed factor in (2.9) as 

,— ind 



(2.10) 



T v2vr 



k{x,x')d6 



^—inO 



-in6' 



T v27T 

e -in(6-6') 

T 



/c(# — , r, z, r , z) dO 



2tt 



k(e-e',r,z,r',z')d9 



-ind' 



k n (r,z,r',z'). 



Inserting (2.10) into (2.9) and executing the integration of 9' over T, we find that (2.1) is 
equivalent to the sequence of equations 

(2.11) a n (r, z) + V2ir k n (r,z,r',z')a n (r',z')r' dl(r',z') = f n (r,z), n£Z. 
For future reference, we define for n £ Z the boundary integral operators K, n via 

(2.12) [Kn a n ] (r, z) = V2n / k n (r,z,r',z')a n (r',z')r' dl(r',z'). 



Then equation (2.11) can be written 
(2.13) 



I + JC n )a n = f n , n e Z. 



When each operator J + /C n is continuously invertible, we can write the solution of (2.1) as 

ind 

(2.14) a(r, z,9) = Y] ~^=[{I + ^)~Vn](n *). 

— * - ' Z7T 



2.3. Truncation of the Fourier series. When evaluating the solution operator (2.14) in 
practice, we will choose a truncation parameter Np, and evaluate only the lowest 2Np + 1 
Fourier modes. If Np is chosen so that the given function / is well-represented by its lowest 
2Np + 1 Fourier modes, then in typical environments the solution obtained by truncating 
the sum (2.14) will also be accurate. To substantiate this claim, suppose that e is a given 
tolerance, and that Np has been chosen so that 

(2.15) 

We define an approximate solution via 
(2.16) 



E ^fn\\<e, 



n=-N v 



2vT 



approx 



ind 



2ir 



n=-N F 

From Parseval's identity, we then find that the error in the solution satisfies 



a - a. 



approx | 



J2 ii(/+/c„rVnii 2 < £ iKi+^rYii/" 



|n|>JV P 



|n|>7V P 



< (max ||(J + /C n )- 1 || 2 ) Yl ll/n|| 2 < (max ||(I + /C B )- 1 || 2> ) e\ 

It is typically the case that the kernel k(x, x') has sufficient smoothness such that the Fourier 
modes k n (r, z, r' , z') decay as n — > oo. Then \\IC n \ \ — > as n — > oo and ||(J + /C n ) _1 || — > 1. 
Thus, an accurate approximation of / leads to an approximation in a that is of the same 
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order of accuracy. Figure 6.3 illustrates that when k is the double layer kernel associated 
with the Laplace equation, and 7 is a simple curve, then + /C n ) _1 || — > 1 with rapid 
convergence. 



3. Discretization of BIEs in two dimensions 

The technique in Section [2] reduces the BIE (1.1) defined on an axisymmetric surface 
r = 7 x T contained in M 3 , to a sequence of BIEs defined on the curve 7 contained in M 2 . 
These equations take the form 



(3.1) 



a(x) + V2tt / k n (x, x') a(x') r' dl(x') = f(x), x G 7, 

J-y 



where the kernel k n is defined as in (2.5). In this section, we describe some standard 
techniques for discretizing an equation such as (3.1). For simplicity, we limit attention to 
the case where 7 is a smooth closed curve, but extensions to non-smooth curves can be 
handled by slight variations of the techniques described here, [U E] • 

3.1. Parameterization of the curve. Let 7 be parameteriz ed b y a vector- valued smooth 
function r : [0, T] — > M 2 . The parameterization converts (3.1) to an integral equation 
defined on the interval [0, T]: 

(3.2) a{r{t))+yfa [ k n (r(t), t(s)) <t(t(s)) r'(r(s)) \dr/ds\ ds = f(r(t)), t G [0, T). 

Jo 

To keep our formulas uncluttered, we suppress the parameterization of the curve and the 
dependence on n and introduce a new kernel 

(3.3) K(t,s) = V2TTk n {r(t),T(s))r'(T(s)) \dr/ds\, 
as well as the functions 

<p(t) = tr(r(t)) and jjj{t) = /(r(t)). 
Then techniques for solving 



(3.4) 



<p(t)+ K(t,s)(p(s)ds = ip(t), t€[0,T\, 



where ip is given and 99 is to be determined, will be equally applicable to (3.2). 



3.2. Nystrom method. We will discretize (3.4) via Nystrom discretization on standard 
Gaussian quadrature nodes, see [2]. To this end, we divide the interval SI = [0, T] into a 
disjoint partition of A^p intervals, 

N P 

n=\Jn p , 

where each £l p is a subinterval called a panel. On each panel O p , we place the nodes of a 

standard A^Q-point Gaussian quadrature rule {t\ The idea is now to enforce (3.4) at 

each of the NpNq nodes: 

a(tS p) ) + f K^^) p(s) ds = ^(tS p) ), (i,p) G {1, 2, . . . , N G } x {1, 2, . . . , iV P }. 
Jo 



To obtain a numerical method, suppose that we can construct for p, q £ {1, 2, . . . , A^p} 
and i,j £ {1, 2, Ag} numbers A^f-^ such that 

T N P N G 

(3.5) / K(t?\s)rts)ds*J2J2 A % 9) ^ ) y 

J ° g=l j=l 

Then the Nystrom method is given by solving the linear system 

Np N G 

(3-6) ^ + EE A tf ] <PT = M e {1, 2, . . . , Ag} x {1, 2, ... , iVp}, 

q=l j = l 



where ip\ = ip(t\ ) and <^ approximates <p(ti ). We write (3.6) compactly as 

(I + A) (p = i/> 

where A is a matrix formed by A^p x A^p blocks, each of size Nq x Nq. We let A^ denote 
the block of A representing the interactions between the panels £l p and £l q . 

3.3. Quadrature and interpolation. We need to determine the numbers A^f- such that 



(3.5) holds. The detailed construction is given in Section 3.4, and utilizes some well-known 



techniques of quadrature and interpolation, which we review in this section. 

3.3.1. Standard Gaussian quadratures. Given an interval [0, h] and a positive integer Ag, 
the Ac-point standard Gaussian quadrature rule consists of a set of A^g nodes {tj}^^ C 



[0, h], and Ag weights {wj}^^ such that 



o 



h N G 

g(s)ds = J2 w j9{tj) 

3=1 



whenever g is a polynomial of degree at most 2A r G — 1, and such that 

,h X G 



[ g(s)ds = Y / w j g(t j ) + 0(h 2N *), 



.]-- 

whenever g is a function with 2A^g continuous derivatives, see pp. 

3.3.2. Quadrature rules for singular functions. Now suppose that given an interval [0, h] 
and a point t £ [—h, 2h], we seek to integrate over [0, h] functions g that take the form 

(3.7) g(s) = log |a-t| + 02(a), 

where cfti and <p2 are polynomials of degree at most 2A^g — 1. Standard Gaussian quadrature 



would be highly inaccurate if applied to integrate (3.7). Rather, seek a A^ G -node quadrature 
that will evaluate 



(3.8) ( g(s)ds 

Jo 

exactly. Techniques for constructing such generalized quadratures are readily available in 
the literature, see for example (15]. These quadratures will be of degree 2A 7 g — 1, just as 
with standard Gaussian quadratures and exhibit comparable accuracy, although in general 
Nq > Ag- The generalized quadratures used in this paper were determined using the 
techniques of [IS], and can be found in the appendix. 

We observe that the quadrature nodes constructed by such methods are typically different 
from the nodes of the standard Gaussian quadrature. This complicates the the construction 



of the matrix A, as described in Section 3.4 



3.3.3. Lagrange interpolation. Let {tj } denote the nodes of a iVc-point Gaussian quad- 
rature rule on the interval [0, h]. If the values of a polynomial g of degree at most Nq — 1 
are specified at these nodes, the entire polynomial g can be recovered via the formula 

N G 
3=1 

where the functions Lj are the Lagrange interpolating polynomials 

s - U 



n 



If g is a smooth function with Nq continuous derivatives that is not a polynomial, then the 
Lagrange interpolant provides an approximation to g satisfying 

N G 

g{s)-Y J L j {^)9{t j ) <Ch N e, 



where 



C=[ sup \g^\s)\ )/N G \. 
\se[o,h] J 



3.4. Constructing the matrix A. Using the t ools reviewed in Section 3.3 



in position to construct numbers Afl- such that (|3.5l) holds. We first note that in forming 



we are now 



i,3 



block Ab'ti of A, we need to find numbers A^f- q ^ such that 



(3.9) f K{t { f\s)^{s)ds^Y.A\ 

J fig j — l 



(P>?) 

3 



<P(t 



1, 2, N G . 



When Q p and Sl g are well separated, the integrand in (3.9) is smooth, and our task is easily 
solved using standard Gaussian quadrature (as described in Section 3.3.1): 

N G 



/ K(t^\s) (f(s) ds *^w j K$\t®)v(t®). 
It directly follows that the ij entry of the block A^ p ' q ^ takes the form 



A 



Wj K(tf\tf). 



(3-10) 

Complications arise when we seek to form a diagonal block A^' p \ or even a block that 
is adjacent to a diagonal block. The difficulty is that the kernel K(t, s) has a singularity as 
s — )■ t. To be precise, for any fixed t, there exist smooth functions ut and vt such that 

K (t, s) = log \t - s\ u t (s) + v t (s). 

(v) \ I 

We see that when t\ is a point in £l q the integrand in (3.9) becomes singular. When 

(p) 

if is a point in a panel neighboring Q q , the problem is less severe, but Gaussian quadra- 
ture would still be inaccur ate. T o maintain full accuracy, we use the modified quadrature 
rules described in Section 3.3. 2[ For every node t^ £ Q p , we construct a quadrature 
{^.tS'^S such that 

(3.11) 



N G 

f K{t ( f\s)v{s)ds^w 
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In order to have a quadrature evaluated at the Gaussian nodes G tt q , we next use 

Lagrange interpolation as described in Section 3.3.3 With {Lj}^^ denoting the Lagrange 
interpolants of order Nq — 1 defined on Q q , we have 



(3.12) 



N G 



t G Q q . 



Inserting (3.12) into (3.11), we find that 



, N a N G 

/ K^\s)^s)ds^Y.^f )K ^ ) ^f ) ) E L ?$f)rt?) 

JSlq £=1 j = l 



We now find that the block A^ p,q ^ of A has entries 
(3.13) ^ = ^*^K{t<*\^)L®$ 



i,j€ {1,2,..., N G }. 



We observe that the formula (3.13) is quite expensive to evaluate; in addition to the sum- 

(p) 

mation, it requires the construction of a quadrature rule for each point Q and evaluation 

of Lagrange interpolants. Fortunately, this process must be executed for at most three 

blocks in each row of blocks of A. 



4. A GENERAL ALGORITHM 



4.1. Summary. At this point, we have shown how to convert a BIE defined on an ax- 
isymmetric surface in M 3 to a sequence of equations defined on a curve in R 2 (Section [2]), 
and then how to discretize each of these reduced equations (Section [3]). Putting everything 
together, we obtain the following algorithm for solving ( |1.1[ ): 

(1) Given the right hand side /, and a computational tolerance e, determine a trunca- 



tion parameter Np such that (2.15) holds. 



(2) Form for 



n 



-Np, —Np + 1, — Np + 2, . . . , Np the matrix A n discretizing the 
equation (2.13) encapsulating the n'th Fourier mode. The matrix is formed via 



Nystrom discretization as described in Section [3] with the discretization parameters 
iVp and Nq chosen to meet the computational tolerance e. 



(3) Evaluate via the FFT the terms {fn}nZ-N m ^ ne Fourier representation of / (as 



defined by (2.3)), and solve for n - 
(I + A n ) a n = f n for a n . Construct a, 
FFT. 



Np, -Np + l, -N F + 2, 



approx 



, A^f the equation 
using formula (2.16) evaluated via the 



The construction of the matrices A n in Step 2 can be accelerated using the FFT (as 



described in Section 4.2), but even with such acceleration, it is typically by a wide margin 



the most expensive part of the algorithm. However, this step needs to be performed only 
once for any given geometry, and given discretization parameters Np, Np, and Nq. The 



method therefore becomes particularly efficient when ( 1.1 ) needs to be solved for a sequence 



of right-hand sides. In this case, it may be worth the cost to pre-compute the inverse of 
each matrix I + A n . 
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4.2. Techniques for forming the matrices. We need to construct for each Fourier mode 
n, a matrix A n consisting of Np x Np blocks An' 9 \ each of size Nq x Nq. Constructing 
an off-diagonal block An ,q ^ when Q p and O g are not directly adjacent is relatively straight- 



forward. For any pair of nodes t^ £ fi p and t^ 



6 Q q , we need to construct the numbers, 



(p>?) 



cf. Q and (|330|), 

for ra = —Np, —Np + 1, . . . , A^p, where r is a parameterization of 7 (see Section 3.1) and 
the kernel k n is defined by (2.5). Fortunately, we do not need to explicitly evaluate the 



' 2-7T Wj k r . 



{r{&\T(tf))r\r{tf))\dr{t^)/ds 



integrals in (2.5) since all the 2Np + 1 numbers can be evaluated by a single application of 



the FFT to the function 
(4.2) 



6^k(6, T(t^),T(tf)). 



(ph 



is not close to T{t)j 



the function in (4.2) is smooth, and the trapezoidal 



When r(t 

rule implicit in applying the FFT is highly accurate. 

Evaluating the blocks on the diagonal, or directly adjacent to the diagonal is somewhat 



more involved. The matrix entries are now given by the formula, cf. (3.3) and (3.13), 



(4.3) 



A 



(p,q) 

k;i,j 



"'a 



„fort) 1, 



(T(t 



(ph 



(?) r forth 



where r and k n are as in (4.1 ). To further complicate things, the points r(t^) and r(t 



forth 



are now in close proximity to each other, and so the functions 



(4.4) 



r(t^),r(i^')) 



tfp,qh 



have a sharp peak around the point 9 = 0. They are typically still easy to integrate away 



from the origin, so the integrals in (2.5) can for a general kernel be evaluated relatively 



efficiently using quadratures that are adaptively refined near the origin. 

Even with the accelerations described in this section, the cost of forming the matrices A n 
tends to dominate the computation whenever the kernels k n must be evaluated via formula 



(2.5). In particular environments, it is possible to side-step this problem by evaluating the 



integral in (2.5) analytically. That this can be done for the single and double layer kernels 



associated with Laplace's equation is demonstrated in Section [5] 

4.3. Computational costs. The asymptotic cost of the algorithm described in Section 
4.1 has three components: (a) the cost of forming the matrices {A n }^_ N , (b) the cost of 
transforming functions from physical space to Fourier space and back, and (c) the cost of 
solving the linear systems (I + A n ) a n = f n . In this section, we investigate the asymptotic 
cost of these steps. We consider a situation where Np Fourier modes need to be resolved, 
and where Np x Nq nodes are used to discretize the curve 7. 

(a) Cost of forming the linear systems: Suppose first that we have an analytic formula for 



each kernel k n . (As we do, e.g., when the original BIE (1.1) involves either the single or 
the double layer kernel associated with Laplace's equation, see Section [5]) Then the cost 
T ma t °f forming the matrices satisfies 

+ Nq Np 



mat 



AfJ Nq Np 

y v ' 

cost from kernel evaluations 



cost from composite quadrature 



When the kernels have to be evaluated numerically via formula (2.5), the cost of forming 



the matrices is still moderate. In the rare situations where the kernel is smooth, standard 
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Gaussian quadrature can be used everywhere and the FFT acceleration described in Section 



4.2 can be used for all entries. In this situation, 

T mat ~N*NlN F log(JV F ). 

In the more typical situation where each kernel k n involves an integrable singularity at the 
diagonal, the FFT acceleration can still be used to rapidly evaluate all entries well-removed 
from the diagonal. However, entries close to the diagonal must be formed via the composite 
quadrature rule combined with numerical evaluation of k n via an adaptive quadrature. In 
this situation, 

T mat ~ N$ N G Np log(iV F ) + N P N G N$. 

(b) Cost of Fourier transforms: The boundary data defined on the surface must be con- 
verted into the Fourier domain. This is executed via the FFT at a cost Tgt satisfying 

(4.5) T m ~NpN G Np log(jV P ). 



We observe that the constant of proportionality in (4.5) is very small, and the cost of this 



step is typically negligible compared to the costs of the other steps. 

(c) Cost of linear solves: Using standard Gaussian elimination, the cost T so i ve of solving 
Ap linear systems (J + A n ) a n = f n , each of size A^p Nq x A^p A^g, satisfies 

r so l ve ~ Ap Nq Np. 

In situations where the equations need to be solved for multiple right hand sides, it pays off 
to first compute the inverses (J + A n ) _1 , and then simply apply these to each right hand 
side (or, alternatively, to form the LU factorizations, and then perform triangular solves). 
The cost T; nv of computing the inverses, and the cost T app iy of applying them then satisfy 

Ti m ~ Np Nq Np , 

tappiy ~ ATp Nq Np . 



We make some practical observations: 

• The cost of forming the matrices by far dominates the other costs unless the kernel 
is either smooth, or analytic formulas for k n are available. 

• The scheme is highly efficient in situations where the same equation needs to be 
solved for a sequence of different right hand sides. Given an additional right hand 
side, the added cost T so i ve is given by 

Toive ~ Ap Nq Np log(AT F ) + Np 1 N G N f , 

with a very small constant of proportionality. We note that this cost remains small 
even if an analytic formula for k n is not available. 

• The system matrices / + A n often have internal structure that allow them to be 
inverted using "fast methods" such as, e.g., those in |17j . The cost of inversion and 
application can then be accelerated to near optimal complexity. 



5. Simplifications for the double layer kernels associated with Laplace's 

equation 

5.1. The double layer kernels of Laplace's equation. Let D C M 3 be a bounded 
domain whose boundary is given by a smooth surface r, let E = D c denote the domain 
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exterior to D, and let n and be the outward unit normal to D. Consider the interior and 
exterior Dirichlet problems of potential theory [IT] , 

(5.1) An = in D, u = f on T, (interior Dirichlet problem) 

(5.2) An = in E, u = f on T. (exterior Dirichlet problem) 



The solutions to (|5.1|) and (|5.2|) can be written in the respective forms 
u(x) 



n(x') ■ (x — x') „ 
v ' y ig-MaO dA(x'), xeD, 



r Att\x — x'\ 



u(x) 



n{x 



(x - x') 



+ 



1 



a(x')dA(x'), x £ E, xq £ D, 



4ir\x — x'\ 3 ' 4ir\x — Xq\ 

where a is a boundary charge distribution that can be determined using the boundary 
conditions. The resulting equations are 



(5.3) 
(5.4) 



<j(x) + 



n(x') • (x — x') 



Mx) + 



Ait\x — x'\ 3 
./ 



nix 



x — x 



a(x')dA(x') = f(x), x£T, 
1 



4tt\x — x'\ 3 



+ 



a(x')dA(x') = f(x), x£T. 



4ir\x — xq\ 

Remark 5.1. There are other integral formulations for the solution to Laplace's equation. 
The double layer formulation presented here is a good choice in that it provides an integral 
operator that leads to well conditioned linear systems. However, the methodology of this 
paper is equally applicable to single-layer formulations that lead to first kind Fredholm 
BIEs. 

5.2. Separation of variables. Using the procedure given in Section [2j if V = 7 x T, then 
(5.1) and (5.2) can be recast as a series of BIEs defined along 7. We express n in cylindrical 



coordinates as 
Further, 



n(x') = (n r > cos 6\n r i s\n9',n z i). 



\x — x 



l\2 



[r cos 

„2 1 / „/\2 



r' cos#') 2 + (rsin( 



r sin6*') 2 + (z 



J\2 



r 2 + (r') 2 — 2rr'(sin 9 sin 0' + cos 6 cos 9') + (z — z') 2 
r 2 + (r') 2 - 2rr' cos((9 - 9') + (z - zf 



and 



n(x') ■ (x — x') = (n r i cos 9' , n r i sin 9' , n z i) ■ (r cos 9 — r cos 9' , r sm.9 — r sin 9' , z — z) 
= n r /r(sin 9 sin 9' + cos 9 cos 9') — n r /r' + n z i (z — z') 
= n r i(r cos{9 — 9') — r) + n z i(z — z). 
Then for a point x' £ T, the kernel of the internal Dirichlet problem can be expanded as 

47T \x — a;' 



where 



z,r ,z 



1 



2tt ^ 



V32tt 3 Jt 



-inO 



n r i (r cos 9 — r') + n z > [z — z') 
(r 2 + (r') 2 - 2rr' cos 9 + (z - z') 2 ) 3 / 2 



d9. 



Similarly, the kernel of the external Dirichlet problem can be written as 



nix 



(x - x') 



4tt\x — x'\ 3 



+ 



1 



4ir\x — Xq\ 



2vr 4^ 



r,z,r ,z 
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with 



4 e) (wV) 



i 



V327T 3 Jf 



—ind 



+ 



n T i (r cos 9 — r') + n z > {z — z') 
(r 2 + (r') 2 - 2rr' cos 9 + (z - z') 2 ) 3 /'- 



dO, 



+ 



(r 2 +Vq — 2ttq cos 9 + (z — zq) 2 ) 1 / 2 



where xq has been written in cylindrical coordinates as (ro cos(#o)) r o sin(#o), zq). With the 
expansions of the kernels available, the procedure described in Section [4] can be used to 
solve (5.3) and (5.4) by solving 



(5.5) 

and 
(5.6) 



a n (r, z) + V2tt / d^(r,r',z,z')a n (r',z')r' dl(r',z') = f n (r,z) 



a n (r,z) + V2^ / d ( n e \r,r',z,z')a n (r',z')r'dl(r',z') = f n (r,z), 



respectively for n - 
log-singularity when both r' = r and z' 



Np, —Np + 1, . . . , Np. Note that the kernels d$ and dffi contain a 



z. 



in cylindrical coordinates, 



Equivalently, (5.5) and (5.6) can be arrived at by considering Laplace's equation written 

0. 



a 2 - 



i d 2 



"u 1 du 
d 2 r r dr r 2 d 2 9 



+ 



2 



u 



d 2 2 



Taking the Fourier transform of u with respect to theta 9 gives 

52„, 1 a„, ^2 a2„, n2„ 



where e n = e n (9) = e m& /-\/2n and u 
with this sequence of PDEs. 



d u n 1 du n n 2 d u n d u n 



r 2 d 2 9 



9 2 z 



0, n G 



e„« n . Then (5.5) and (5.6) are now associated 



5.3. Evaluation of kernels. The values of dn and for n = —Np, —Np+1, . . . , A^p need 

to be com pute d efficiently and with high accuracy to construct the Nystrom discretization of 

I I I (i) ( e ) 

(5.5) and (5.6). Note that the integrands of dn and d n are real valued and even functions 

on the interval [— tv,tt]. Therefore, d$ can be written as 



(5.7) d<Z\r,z,r',z') 



1 



n r i (r cos t — r') + n z i {z — z') 
(r 2 + (r') 2 - 2rr' cos t + (z- z') 2 ) 3 / 2 



cos(ni) dt. 



(e) 

Note that d n can be written in a similar form. 

This integrand is oscillatory and increasingly peaked at the origin as both r' — > r and 
z' —> z. As long as r' and r as well as z 1 and z are well separated, the integrand does not 
experience peaks near the origin, and as mentioned before, the FFT provides a fast and 
accurate way for calculating dn and dn\ 

In regimes where the integrand is peaked, the FFT no longer provides a means of evalu- 
ating dn and dn" 1 with the desired accuracy. One possible solution to this issue is applying 
adaptive quadrature to fully resolve the peak. However, this must be done for each value 
of n required and becomes prohibitively expensive if Np is large. 
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Fortunately, an analytical solution to (5.7) exists. As noted in [6], the single-layer kernel 
can be expanded with respect to the azimuthal variable as 



s(x, x') 



A-k\x - x'\ 4ti-2 



nez 



where Q n -\/2 is the half-integer degree Legendre function of the second kind and 

_ r 2 + {r') 2 + (z - z') 2 
~^ 2rr' 

In light of this expansion, single-layer kernel can similarly be written as 
s(x, X ) 



4vr(r 2 + (r') 2 - 2rr> cos(<9 - 9') + (z - z') 2 ) 1 / 2 
-^^e^') Sn (r,,,r',z') 

V Z7T — 



where 



s n (r,z,r',z') 



1 



cos(nt) 



\/32^3 Jj (r 2 + (r') 2 - 2rr< cos(t) + (z - z') 2 ) 1 / 2 
1 / cos(ni) 



V 87r 3 rr' Jt \/8(x — cos(t)) 

=== Qn-l/2(x) 
V ovr^rr' 



To find an analytical form for (5.7), first note that in cylindrical coordinates the double- 



layer kernel can be written in terms of the single-layer kernel, 

n(x') ■ (x — x') n r / (r cos(6> — 9') — r') + n z > (z — z') 



4tt\x - x'\ 3 4vr(r 2 + (r') 2 - 2rr' cos{9 - 6') + (z - z') 2 ) 3/2 



1 

47T 

+ n z 



d_ f 1 

d? V (r 2 + (r') 2 - 2rr' cos(# - 9') + (z - z') 2 ) 1 / 2 

d_ ( 1 

dz 1 V (r 2 + (r') 2 - 2rr' cos(9 - 9') + (z - z') 2 ) 1 / 2 



+ 



The coefficients of the Fourier series expansion of the double-layer kernel are then given by 

(i) 

dn , which can be written using the previous equation as 

d ( cos(nt) 



d$(r,z,r',z') =n r < 



+ n 



T dr' V (32vr 3 (r 2 + (r') 2 - 2rr' cos(i) + (z - z') 2 )) 1 / 2 
d ( cos(ni) 



T 



dz' V (32^ 3 (r 2 + (r') 2 - 2rr' cos(t) + {z- z') 2 )) 1 / 2 



dt+ 
dt 



d 

• 

dr 
1 



fdQ n - 1/2 (x)d X Qn-1/2(X)\ dQ n _ 1/2 ( X )d X 
n r i I -tt^ —. I + n z i 



dx dr' 



dx dz' 



15 



To utilize this form of dn , set /x = J — and note that 



dx {r') 2 — r 2 — (z — z') 2 

dr' 2r(r / ) 2 

dx z' — z 

dz' rr' ' 

Q-i/ 2 (x)=^W, 

Qi/ 2 (X) = XV-Kiv) - ^2{ X + l)E(n), 

2-n-l/ 2 (x) = Qn-1/2(X), 

Sn-1/2(X) = 4 2n - 1 ^Q"~3/ 2 (X) - 2n _ 1 Qn-5/2(x)> 

<9Qn-i/ 2 (x) 2n- 1 , , 
~1?X ~~ 2( X 2 - 1) ^ xyn ^ 1 / 2 ~ y "~3/2J ' 



where -ff and E are the complete elliptic integrals of the first and second kinds, respectively. 
The first two relations follow immediately from the definition of x an d the relations for the 
Legendre functions of the second kind can be found in pQ. With these relations in hand, 

(i) 

the calculation of d n for n = —Np, —Np + 1, . . . , Np can be done accurately and efficiently 
when r 1 and r as well as z' and z are in close proximity. The calculation of dn" 1 can be done 
analogously. 

Remark 5.2. Note that the forward recursion relation for the Legendre functions Q n _i/2(x) 
is unstable when % > 1. In practice, the instability is mild when x is near 1 and can still 
be employed to accurately compute values in this regime. Additionally, if stability be- 
comes and issue, Miller's algorithm [H] can be used to calculate the values of the Legendre 
functions using the backwards recursion relation, which is stable for x > !■ 



6. Numerical results 
This section describes several numerical experiments performed to assess the efficiency 



and accuracy of the the numerical scheme outlined in Section 4.1 All experiments were 



executed for the double layer kernels associated with Laplace's equation, calculated using 



the recursion relation described in Section 5.3. Note that the kernels in this case give us the 



property that A- n = A n , and so we need only to invert A^f + 1 matrices. Further, the FFT 
used here is complex-valued, and a real-valued FFT would yield a significant decrease in 
computation time. The geometries investigated are described in Figure [67T| The generating 
curves were parameterized by arc length, and split into Ap equisized panels. A 10-point 
Gaussian quadrature has been used along each panel, with the modified quadratures used 
to handle the integrable singularities in the kernel. These quadratures are listed in the 
appendix. The algorithm was implemented in MatLab and the experiments were run on a 
2.66GHz Intel Quad Core with 6Gb of RAM. All timings were averaged over 10 runs. 



6.1. Computational costs. Using the domain in Figure [6?l| a) and the interior Dirichlet 
problem, timing results are given in Table 1. The reported results include: 
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(a) 





(b) 




(c) 

Figure 6.1. Domains used in numerical examples. All items are rotated 
about the vertical axis, (a) A sphere, (b) A wavy block, (c) A starfish 
torus. 



Np 
N F 

?setup 
^mat 

T- 

^apply 



the number of panels used to discretize the contour 

the Fourier truncation parameter (we keep 2Np + 1 modes) 

time to setup the discrete system, excluding construction of the linear systems 

time to construct the linear systems (utilizing the recursion relation) 

time to invert the linear systems 

time to Fourier transform the right hand side and the solution 
time to apply the inverse to the right hand side 



The most expensive component of the calculation is the construction of the linear systems. 
This is primarily a result of the cost of evaluating the kernel and applying the modified 
quadrature rules. Table 2 compares the use of the recursion relation in evaluating the kernel 
when it is near-singular to using an adaptive Gaussian quadrature. The efficiency of the 



recursion relation is clearly evident in this case. Figure 6.2 plots the time to construct the 
linear systems as the number of panels and as the number of Fourier modes increases. The 
time to construct the systems as the number of panels increases grows as O(NpNq) as 
expected, and the timings go as O(Np) as the number of Fourier modes increase. Note that 
it is rather inexpensive to increase the number of Fourier modes used, making the scheme 
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particularly attractive in environments where a large number is required. As an artifact 
of the small lengths of the vectors being Fourier transformed, the asymptotic behavior of 
the timings as the number of Fourier modes increases is not evident until we reach a large 
number of modes. 

We observe that the largest problem reported in Table 1 involves 320 000 degrees of 
freedom. The method requires 2.2 minutes of pre-computation for this example, and is 
then capable of computing a solution u from a given data function / in 0.46 seconds. 



7Vp 


2Np + 1 


-^setup 


^mat 




Tm 


^apply 


5 


25 


9.10e-02 


4.31e-01 


2.48e-03 


1.17e-03 


3.82e-04 


10 


25 


1.10e-01 


l.lye+UU 


n o 1 „ no 


1 flQn no 

l.yoe-Uo 


r r)-i r\i 

5.zle-04 


20 


25 


1.22e-01 


3.42e+00 


4.47e-02 


3.42e-03 


1.20e-03 


40 


25 


1.68e-01 


1.15e+01 


2.90e-01 


6.66e-03 


4.94e-03 


80 


25 


2.60e-01 


4.90e+01 


2.12e+00 


1.32e-02 


1.94e-02 


5 


50 


9.63e-02 


4.59e-01 


4.79e-03 


1.63e-03 


6.85e-04 


10 


50 


l.lbe-01 


1 0^7„ i f\C\ 

1.27e+00 


1.7ze-l)z 


o oi „ no 

2.81e-03 


9.71e-04 


20 


50 


1.54e-01 


3.99e+00 


8.89e-02 


5.47e-03 


2.86e-03 


40 


50 


2.32e-01 


1.34e+01 


5.80e-01 


1.08e-02 


1.00e-02 


80 


50 


4.04e-01 


5.05e+01 


4.27e+00 


2.13e-02 


3.89e-02 


5 


100 


1.13e-01 


4.72e-01 


9.20e-03 


2.74e-03 


1.29e-03 


1 n 

1U 


1UU 


i.^iye-ui 


i ^/ic._i_nn 
i.o^ie-i-uu 


o.o4e-uz 


^.yoe-uo 


i .yye-uo 


20 


100 


2.20e-01 


4.19e+00 


1.75e-01 


9.80e-03 


6.20e-03 


40 


100 


3.75e-01 


1.47e+01 


1.16e+00 


2.06e-02 


2.01e-02 


80 


100 


6.53e-01 


5.76e+01 


8.34e+00 


4.13e-02 


7.99e-02 


5 


200 


1.48e-01 


5.46e-01 


1.92e-02 


5.14e-03 


2.68e-03 


10 


200 


2.24e-01 


1.53e+00 


6.89e-02 


9.53e-03 


4.48e-03 


20 


200 


3.59e-01 


4.90e+00 


3.44e-01 


1.88e-02 


1.25e-02 


40 


200 


6.45e-01 


1.63e+01 


2.29e+00 


3.74e-02 


4.13e-02 


80 


200 


1.22e+00 


6.85e+01 


1.68e+01 


7.67e-02 


1.58e-01 


5 


400 


2.11e-01 


8.49e-01 


3.81e-02 


1.01e-02 


5.46e-03 


10 


400 


3.45e-01 


2.15e+00 


1.29e-01 


1.92e-02 


9.59e-03 


20 


400 


5.83e-01 


6.47e+00 


6.85e-01 


3.81e-02 


2.54e-02 


40 


400 


1.09e+00 


2.22e+01 


4.49e+00 


7.60e-02 


7.94e-02 


80 


400 


2.27e+00 


9.33e+01 


3.27e+01 


1.51e-01 


3.06e-01 



Table 1. Timing results in seconds performed for the domain given in 
Figure |6.1|^a) for the interior Dirichlet problem. 



6.2. Accuracy and conditioning of discretization. The accuracy of the discretization 
has been tested using the interior and exterior Dirichlet problems on the domains given in 



Figure 6.1 Exact solutions were generated by placing a few random point charges outside 
of the domain where the solution was calculated. The solution was evaluated at points 
defined on a sphere encompassing (or interior to) the boundary. The errors reported in 
Tables 3, 4, and 5 are relative errors measured in the l°° -norm, 1 1 u e — ^||oo/||^||oo? where u 
is the exact potential and u e is the potential obtained from the numerical solution. 

In all cases, 10 digits of accuracy has been obtained from a discretization involving a 
relatively small number of panels, due to the rapid convergence of the Gaussian quadrature. 
This is especially advantageous, as the most expensive component of the algorithm is the 
construction of the linear systems, the majority of the cost being directly related to the 
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2N F + 1 


Composite Quadrature 


Recursion Relation 


25 


1.9 


0.43 


50 


3.1 


0.46 


100 


6.6 


0.47 


200 


18.9 


0.55 



Table 2. Timing comparison in seconds for constructing the matrices 
(I + A n ) using composite Gaussian quadrature and the recursion relation 



described in Section 5.3 to evaluate k n for diagonal and near diagonal blocks. 
The FFT is used to evaluate k n at all other entries. 27Vp + 1 is the total 
number of Fourier modes used. 5 panels were used to discretize the bound- 
ary. 



mat 




mat 



10 



10 



10 



10 



♦ » ▼ * 




10 



10" 



1 



• 5 

• 10 
20 

» 40 

80 



10 



N P 2N F 
(a) (b) 

Figure 6.2. Timings for the construction of the linear systems, (a) Scaling 
as the number of panels increases. The solid line corresponds to growth with 
the number of panels squared. The legend refers to the number of Fourier 
modes used, (b) Scaling as the number of Fourier modes increases. The 
solid line corresponds to linear growth. The legend refers to the number of 
panels used. 



number of panels used. Further, the number of Fourier modes required to obtain 10 digits of 
accuracy is on the order of 100 modes. Although not investigated here, the discretization 
technique naturally lends itself to nonuniform refinement of the surface, allowing one to 
resolve features of the surface that require finer resolution. 

The number of correct digits obtained as the number of panels and number of Fourier 
modes increases eventually stalls. This is a result of a loss of precision in determining the 
kernels, as well as cancelation errors incurred when evaluating interactions between nearby 
points. This is especially prominent with the use of Gaussian quadratures, as points cluster 
near the ends of the panels. If more digits are required, high precision arithmetic can be 
employed in the setup phase of the algorithm. 

Figure 6.3 shows the singular values as well as the condition numbers of an 80 panel 
discretization for the Np = —200, . . . , 200 Fourier modes used in the discretization of the 
interior Dirichlet problem, on the domain shown in Figure [(Tl^ a). The integral equations 
of this paper are second kind Fredholm equations, and generally lead to well-conditioned 
systems. As seen in Figure 6.3, this hold true for the discretization presented in this paper. 
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N P 


2N F + 1 




25 


50 


100 


200 


400 


5 


1.93869e-04 


4.10935e-07 


5.37883e-08 


5.37880e-08 


5.37880e-08 


10 


1.93869e-04 


4.10513e-07 


3.27169e-12 


6.72270e-13 


6.72270e-13 


20 


1.93869e-04 


4.10513e-07 


3.30601e-12 


1.66132e-13 


1.66132e-13 


40 


1.93869e-04 


4.10513e-07 


3.23162e-12 


8.28568e-14 


8.28568e-14 


80 


1.93869e-04 


4.10512e-07 


2.92918e-12 


2.92091e-13 


2.92091e-13 



Table 3. Error in internal Dirichlet problem solved on domain (a) in Figure 6.1 



N P 


2N F + 1 




25 


50 


100 


200 


400 


5 


9.11452e-04 


9.11464e-04 


9.11464e-04 


9.11464e-04 


9.11464e-04 


10 


4.15377e-05 


4.15416e-05 


4.15416e-05 


4.15416e-05 


4.15416e-05 


20 


6.31923e-07 


1.29234e-07 


1.29235e-07 


1.29235e-07 


1.29235e-07 


40 


7.04741e-07 


3.10049e-ll 


3.08152e-ll 


3.08305e-ll 


3.08359e-ll 


80 


7.04779e-07 


5.62558e-ll 


5.05306e-ll 


5.05257e-ll 


5.05232e-ll 



Table 4. Error in external Dirichlet problem solved on domain (b) in Figure 6.1 



N P 


2N F + 1 




25 


50 


100 


200 


400 


5 


3.80837e-04 


3.83707e-04 


3.83707e-04 


3.83707e-04 


3.83707e-04 


10 


2.41602e-05 


6.81564e-06 


6.81556e-06 


6.81556e-06 


6.81556e-06 


20 


3.03272e-05 


5.98506e-09 


2.53980e-ll 


2.54112e-ll 


2.54118e-ll 


40 


3.03272e-05 


6.01273e-09 


6.95662e-12 


6.94592e-12 


6.94546e-12 


80 


3.03272e-05 


6.01059e-09 


5.25217e-12 


5.26674e-12 


5.26515e-12 



Table 5. Error in external Dirichlet problem solved on domain (c) in Figure 6.1 



7. Generalizations and Conclusions 

This paper describes a numerical technique for computing solutions to boundary inte- 
gral equations defined on axisymmetric surfaces in IR 3 with no assumption on the loads 
being axisymmetric. The technique is introduced as a generic method with only very mild 
conditions imposed on the kernel; specifically, we assume that the kernel has an integrable 



singularity at the diagonal, and that it is rotationally symmetric (in the sense that (2.2) 
holds). The technique described improves upon previous work in two regards: 

(1) A highly accurate quadrature scheme for kernels with integrable singularities is 
introduced. Numerical experiments indicate that solutions with a relative accuracy 
of 10~ 10 or better can easily be constructed. 



(2) A rapid technique for numerically constructing the kernel functions k n in (1.3) 
is introduced. It works when k is either the single or the double layer potential 
associated with Laplace's equation. The technique is a hybrid scheme that relies on 
the FFT when possible, and uses recursion relations for Legendre functions when 
not. The resulting scheme is fast enough that a problem involving 320 000 degrees of 
freedom can be solved in 2.2 minutes on a standard desktop PC. Once one problem 
has been solved, additional right hand sides can be processed in 0.46 seconds. 
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Figure 6.3. Maximum and minimum singular values for the matrices re- 
sulting from an 80 panel discretization of a sphere using 400 Fourier modes, 
where n is the the matrix associated with the n th Fourier mode. 

Some possible extensions of this work include: (1) Acceleration of the solution of the 
linear systems using fast methods, (2) extension of recursion relation to the Helmholtz 
equation, and (3) extension of the algorithm to problems involving multiple bodies whose 
axes of symmetry are not necessarily aligned. 
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APPENDIX OF QUADRATURE NODES AND WEIGHTS 



10 Point Gauss-Legendre Rule for 
integrals of the form f(x) dx 



NODES 



-9.739065285171716e-01 
-8.650633666889845e-01 
-6.794095682990244e-01 
-4.333953941292472e-01 
-1.488743389816312e-01 
1.488743389816312e-01 
4.333953941292472e-01 
6.794095682990244e-01 
8.650633666889845e-01 
9.739065285171716e-01 



WEIGHTS 



6.667134430868814e-02 
1.494513491505806e-01 
2.190863625159820e-01 
2.692667193099963e-01 
2.955242247147529e-01 
2.955242247147529e-01 
2.692667193099963e-01 
2.190863625159820e-01 
1.494513491505806e-01 
6.667134430868814e-02 



20 point quadrature rule for integrals 
of the form f f(x) + g(x) log \xi — x\ dx, 
where xi is a Gauss-Legendre node 
NODES n WEIGHTS 



-9.981629455677877e-01 
-9.915520723139890e-01 
-9.832812993252168e-01 
-9.767801773920733e-01 
-9.717169387169078e-01 
-9.510630103726074e-01 
-9.075765988474132e-01 
-8.382582352569804e-01 
-7.408522006801963e-01 
-6.147619568252419e-01 
-4.615244999958006e-01 
-2.849772954295424e-01 
-9. 1 1759346048974 7e-02 
1.119089520342051e-01 
3.148842536644393e-01 
5.075733846631832e-01 
6.797470718157004e-01 
8.218833662202629e-01 
9.258924858821892e-01 
9.857595961761246e-01 



4.550772157144354e-03 
8.062764683328619e-03 
7.845621096866406e-03 
4.375212351185101e-03 
1.021414662954223e-02 
3.157199356768625e-02 
5.592493151946541e-02 
8.310260847601852e-02 
1.118164522164500e-01 
1.401105427713687e-01 
1.657233639623953e-01 
1.863566566231937e-01 
1.999093145144455e-01 
2.046841584582030e-01 
1.995580161940930e-01 
1.841025430283230e-01 
1.586456191174843e-01 
1. 242680229936 124e-01 
8.273794370795576e-02 
3.643931593123844e-02 



20 point quadrature rule for integrals 
of the form f(x) + g(x) log \x2 — x\ dx, 
where X2 is a Gauss-Legendre node 
NODES I WEIGHTS 



20 point quadrature rule for integrals 
of the form f(x) + g(x) log \x 3 — x\ dx, 
where x$ is a Gauss-Legendre node 
NODES I WEIGHTS 



-9.954896691005256e-01 
-9.775532683688947e-01 
-9.500346715183706e-01 
-9.192373372373420e-01 
-8.916563772395616e-01 
-8.727728136507039e-01 
-8.607963163061316e-01 
-8.201318720954396e-01 
-7.394732321355052e-01 
-6.204853512352519e-01 
-4.667290485167077e-01 
-2.840823320902124e-01 
-8.079364608026202e-02 
1.328455136645940e-01 
3.451233500669768e-01 
5.437321547508867e-01 
7.167077216635750e-01 
8.534299232009863e-01 
9.458275339169444e-01 
9.912353127269481e-01 



1.141744473788874e-02 
2.368593568061651e-02 
3.027205199814611e-02 
3.021809354380292e-02 
2.397183723558556e-02 
1.253574079839078e-02 
2.070840476545303e-02 
6.080709508468810e-02 
1.002402801599464e-01 
1.371499151597280e-01 
1.693838059093582e-01 
1.945292086962893e-01 
2.103223087093422e-01 
2. 14990092844 7852e-01 
2.074984762344433e-01 
1.877085225595498e-01 
1.564543949958065e-01 
1.156104890379952e-01 
6.859369195724087e-02 
2.390220989094312e-02 



-9.930122613589740e-01 
-9.643941806993207e-01 
-9.175869559770760e-01 
-8.596474181980754e-01 
-7.990442708271941e-01 
-7.443700671611690e-01 
-7.031684479828371e-01 
-6.811221147275545e-01 
-6.579449960254029e-01 
-5.949471688137100e-01 
-4.893032793226841e-01 
-3.441659232382107e-01 
-1.665388322404095e-01 
3.344207582228461e-02 
2.434356263087524e-01 
4.498696863725133e-01 
6.389777518528792e-01 
7.978632877793501e-01 
9.155180703268415e-01 
9.837258757826489e-01 



1.779185041193254e-02 
3.870503119897836e-02 
5.371 120494602663e-02 
6.073467932536858e-02 
5.901993373645797e-02 
4.905519963921684e-02 
3.249237036645046e-02 
1.335394660596527e-02 
4.151626407911676e-02 
8.451456165895121e-02 
1.262522607368499e-01 
1.628408264966550e-01 
1.907085686614375e-01 
2.071802230953481e-01 
2.105274833603497e-01 
2.000282912446872e-01 
1.760212445284564e-01 
1.399000904426490e-01 
9.402669072995991e-02 
4.161927873514264e-02 
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20 point quadrature rule for integrals 
of the form J_ 1 f(x) + g(x) log \x± — x\ dx, 
where xa is a Gauss-Legendre node 



20 point quadrature rule for integrals 
of the form f(x) + g(x) log \x$ — x\ dx, 
where x 5 is a Gauss-Legendre node 



NODES 



WEIGHTS 



NODES 



WEIGHTS 



-9.903478871 133073e-01 
-9.504025146897784e-01 
-8.834986023815121e-01 
-7.974523551287549e-01 
-7.022255002503461e-01 
-6.087194789244920e-01 
-5.275278952351541e-01 
-4.677586540799037e-01 
-4.360689210457623e-01 
-4.121945474875853e-01 
-3.49422676691 1471e-01 
-2.425993523586304e-01 
-9.646839923908594e-02 
7.921243716767302e-02 
2.715178194484646e-01 
4.658440358656903e-01 
6.472213975763533e-01 
8.015601619414859e-01 
9.168056007307982e-01 
9.839468743284722e-01 



2.462513260640712e-02 
5.449201732062665e-02 
7.799498604905293e-02 
9.241688894090601e-02 
9.619882322938848e-02 
8.902783806614303e-02 
7.181973054766198e-02 
4.663017060126023e-02 
1.794303974050253e-02 
4.061799823415495e-02 
8.507517518447759e-02 
1.277525783357134e-01 
1.628510773009247e-01 
1.863323765408308e-01 
1.958227701927855e-01 
1.903138548150517e-01 
1.700731513381802e-01 
1.365784674773513e-01 
9.239595239693155e-02 
4.103797108164931e-02 



-9.883561797860961e-01 
-9.398305159297058e-01 
-8.572399919019390e-01 
-7.482086250804679e-01 
-6.228514167093102e-01 
-4.928317114329241e-01 
-3.702771193724617e-01 
-2.666412108172461e-01 
-1.916083010783277e-01 
-1.521937160593461e-01 
-1.233125650067164e-01 
-5.257959675044444e-02 
5.877314311857769e-02 
2.012559739993003e-01 
3.627988191760868e-01 
5.297121321076323e-01 
6.878399330187783e-01 
8.237603202215137e-01 
9.259297297557394e-01 
9.856881498392895e-01 



2.974603958509255e-02 
6.657945456889164e-02 
9.731775484182564e-02 
1.190433988432928e-01 
1.297088242013777e-01 
1.282900896966494e-01 
1.148917968875341e-01 
9.074932908233864e-02 
5.818196361216740e-02 
2.224697059733435e-02 
4.788826761346366e-02 
9.237500180593534e-02 
1.287410543031414e-01 
1.541960911507042e-01 
1.665885274544506e-01 
1.648585116745725e-01 
1.491408089644010e-01 
1.207592726093190e-01 
8.212177982524418e-02 
3.657506268226379e-02 



20 point quadrature rule for integrals 




20 point quadrature rule for integrals 


of the form fix) + g(x) log \xq — x\ dx, 




of the form f(x) + g(x) log \x7 — x\ dx, 


where x& is a Gauss- 


Legendre node 




where x*i is a Gauss- 


Legendre node 


NODES 


WEIGHTS 


NODES 


WEIGHTS 


-9.856881498392895e-01 


3 


657506268226379e- 


02 


-9.839468743284722e-01 


4 


103797108164931e-02 


-9.259297297557394e-01 


8 


212177982524418e- 


02 


-9.168056007307982e-01 


9 


239595239693155e-02 


-8.237603202215137e-01 


1 


207592726093190e- 


01 


-8.015601619414859e-01 


1 


365784674773513e-01 


-6.878399330187783e-01 


1 


491408089644010e- 


01 


-6.472213975763533e-01 


1 


700731513381802e-01 


-5.297121321076323e-01 


1 


648585116745725c- 


01 


-4.658440358656903e-01 


1 


903138548150517e-01 


-3.627988191760868e-01 


1 


665885274544506e- 


01 


-2.715178194484646e-01 


1 


958227701927855e-01 


-2.012559739993003e-01 


1 


541960911507042e- 


01 


-7.921243716767302e-02 


1 


863323765408308e-01 


-5.877314311857769e-02 


1 


287410543031414e- 


01 


9.646839923908594e-02 


1 


628510773009247e-01 


5.257959675044444e-02 


9 


237500180593534e- 


02 


2.425993523586304e-01 


1 


277525783357134e-01 


1.233125650067164e-01 


4 


788826761346366e- 


02 


3.494226766911471e-01 


8 


507517518447759e-02 


1.521937160593461e-01 


2 


224697059733435e- 


02 


4.121945474875853e-01 


4 


061799823415495e-02 


1.916083010783277e-01 


5 


818196361216740e- 


02 


4.360689210457623e-01 


1 


794303974050253e-02 


2.666412108172461e-01 


9 


074932908233864e- 


02 


4.677586540799037e-01 


4 


663017060126023e-02 


3.702771 193724617e-01 


1 


148917968875341e- 


01 


5.275278952351541e-01 


7 


181973054766198e-02 


4.928317114329241e-01 


1 


282900896966494e- 


01 


6.087194789244920e-01 


8 


902783806614303e-02 


6.228514167093102e-01 


1 


297088242013777e- 


01 


7.022255002503461e-01 


9 


619882322938848e-02 


7.482086250804679e-01 


1 


190433988432928e- 


01 


7.974523551287549e-01 


9 


241688894090601e-02 


8.572399919019390e-01 


9 


731775484182564e- 


02 


8.834986023815121e-01 


7 


799498604905293e-02 


9.398305159297058e-01 


6 


657945456889164e- 


02 


9.504025146897784e-01 


5 


449201732062665e-02 


9.883561797860961e-01 


2 


974603958509255e- 


02 


9.903478871 133073e-01 


2 


462513260640712e-02 
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20 point quadrature rule for integrals 




20 point quadrature rule for integrals 


of the form |_j f(x) + g(x) log \x$ — x\ dx, 




of the form f(x) + g(x) log \xg — x\ dx, 


where xg is a Gauss-Legendre node 




where xg is a Gauss-Legendre node 


NODES 


WEIGHTS 


NODES 


WEIGHTS 


-9.837258757826489e-01 


4.161927873514264e- 


02 


-9.912353127269481e-01 


2.390220989094312e-02 


-9.155180703268415e-01 


9.402669072995991e- 


02 


-9.458275339169444e-01 


6.859369195724087e-02 


-7.978632877793501e-01 


1.399000904426490e- 


01 


-8.534299232009863e-01 


1.156104890379952e-01 


-6.389777518528792e-01 


1.760212445284564e- 


01 


-7.167077216635750e-01 


1.564543949958065e-01 


-4.498696863725133e-01 


2.000282912446872e- 


01 


-5.437321547508867e-01 


1.877085225595498e-01 


-2.434356263087524e-01 


2.105274833603497e- 


01 


-3.451233500669768e-01 


2.074984762344433e-01 


-3.344207582228461e-02 


2.071802230953481e- 


01 


-1.328455136645940e-01 


2. 14990092844 7852e-01 


1.665388322404095e-01 


1.907085686614375e- 


01 


8.079364608026202e-02 


2.103223087093422e-01 


3.441659232382107e-01 


1.628408264966550e- 


01 


2.840823320902124e-01 


1.945292086962893e-01 


4.893032793226841e-01 


1.262522607368499e- 


01 


4.667290485167077e-01 


1.693838059093582e-01 


5.949471688137100e-01 


8.451456165895121e- 


02 


6.204853512352519e-01 


1.371499151597280e-01 


6.579449960254029e-01 


4.151626407911676e- 


02 


7.394732321355052e-01 


1.002402801599464e-01 


6.811221147275545e-01 


1.335394660596527e- 


02 


8.201318720954396e-01 


6.080709508468810e-02 


7.031684479828371e-01 


3.249237036645046e- 


02 


8.607963163061316e-01 


2.070840476545303e-02 


7.443700671611690e-01 


4.905519963921684e- 


02 


8.727728136507039e-01 


1.253574079839078e-02 


7.990442708271941e-01 


5.901993373645797e- 


02 


8.916563772395616e-01 


2.397183723558556e-02 


8.596474181980754e-01 


6.073467932536858e- 


02 


9.192373372373420e-01 


3.021809354380292e-02 


9.175869559770760e-01 


5.371 120494602663e- 


02 


9.500346715183706e-01 


3.027205199814611e-02 


9.643941806993207e-01 


3.870503119897836c- 


02 


9.775532683688947e-01 


2.368593568061651e-02 


9.930122613589740e-01 


1.779185041193254e-02 


9.954896691005256e-01 


1.141744473788874e-02 



20 point quadrature rule for integrals 
of the form f_ 1 f(x) + g(x) log \xio — x\ dx, 
where xio is a Gauss-Legendre node 
NODES n WEIGHTS 



-9.857595961761246e-01 
-9.258924858821892e-01 
-8.218833662202629e-01 
-6.797470718157004e-01 
-5.075733846631832e-01 
-3.148842536644393e-01 
-1.119089520342051e-01 
9.117593460489747e-02 
2.849772954295424e-01 
4.615244999958006e-01 
6.147619568252419e-01 
7.408522006801963e-01 
8.382582352569804e-01 
9.075765988474132e-01 
9.510630103726074e-01 
9.717169387169078e-01 
9.767801773920733e-01 
9.832812993252168e-01 
9.915520723139890e-01 
9.981629455677877e-01 



3.643931593123844e-02 
8.273794370795576e-02 
1.242680229936124e-01 
1.586456191174843e-01 
1.841025430283230e-01 
1.995580161940930e-01 
2.046841584582030e-01 
1.999093145144455e-01 
1.863566566231937e-01 
1.657233639623953e-01 
1.401105427713687e-01 
1.118164522164500e-01 
8.310260847601852e-02 
5.592493151946541e-02 
3.157199356768625e-02 
1.021414662954223e-02 
4.375212351185101e-03 
7.845621096866406e-03 
8.062764683328619e-03 
4.550772157144354e-03 



24 point quadrature rule for integrals 
of the form J Q f(x) + g(x) \og(x + x)dx, 
where x > 10 _1 



NODES 



3.916216329415252e-02 
8.135233983530081e-02 
1.123448211344994e-01 
1.595931983965030e-01 
2.085759027831349e-01 
2.426241962027560e-01 
2.886190312538522e-01 
3.469021762354675e-01 
4.072910101569611e-01 
4.664019722595442e-01 
5.182120817844112e-01 
5.501308436771654e-01 
5.970302980854608e-01 
6.548457960388209e-01 
7.119542126106005e-01 
7.607920420946340e-01 
7.953017051155684e-01 
8.303900341517088e-01 
8.612724919009394e-01 
8.954049128027080e-01 
9.315909369155358e-01 
9.621742249068356e-01 
9.843663446380599e-01 
9.970087425823398e-01 



WEIGHTS 



4.880755296918116e-02 
3. 19600278516361 le-02 
3.883416642507362e-02 
5.148898992140820e-02 
4.219328148763533e-02 
3.420686213633789e-02 
5.512488680719239e-02 
6.007112809843418e-02 
6.022350479415180e-02 
5.735022004401478e-02 
4. 1679234171 18068e-02 
3.346089628879600e-02 
5.574716218423796e-02 
5.847838243344473e-02 
5.464156990092474e-02 
4.092186343704961e-02 
3.283728166050225e-02 
3.438233273473095e-02 
3.022585192226418e-02 
3.700769701277380e-02 
3.410213679365162e-02 
2.665791885274193e-02 
1.754420526360429e-02 
7.662283104388867e-03 
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24 point quadrature rule for integrals 
of the form f Q f(x) + g(x) \og(x + x)dx, 
where < x < 10" 1 



NODES 



1.940564616937581e-02 
4.545433992382339e-02 
7.378866604396420e-02 
1.054147718077606e-01 
1.412997888401000e-01 
1.822325567811081e-01 
2.287282121202408e-01 
2.809170925514041e-01 
3.384320962237970e-01 
4.003108031244078e-01 
4.648605571606025e-01 
5.290714994276687e-01 
5.829663557386375e-01 
6.128301889979477e-01 
6.606072156240962e-01 
7.139495966128518e-01 
7.677830914961244e-01 
8.187382423336450e-01 
8.587068551739496e-01 
8.906873285570645e-01 
9.267772492129903e-01 
9.592137652582382e-01 
9.830962712794008e-01 
9.967621546194148e-01 



WEIGHTS 



2.514022176052795e-02 
2.703526530535647e-02 
2.980872487617485e-02 
3.360626237885489e-02 
3.829678083416609e-02 
4.365651045780837e-02 
4.935846322319046e-02 
5.495967924055210e-02 
5.991162198705084e-02 
6.356960862248889e-02 
6.5068685524671 18e-02 
6.219588235225894e-02 
3.889986041695310e-02 
3.573431931940621e-02 
5.296315368353523e-02 
5.369033999927759e-02 
5.340793573367282e-02 
4.704756013998560e-02 
3.276576301747068e-02 
3.449175311880027e-02 
3.560168848238671e-02 
2.857367151127661e-02 
1.894042942442201e-02 
8.291994770212826e-03 



24 point quadrature rule for integrals 
of the form J Q f(x) + g(x) \og(x + x)dx, 
where 1CT 3 < x < 10" 2 



NODES 



7.571097817272427e-03 
1.800655325976786e-02 
3.003901004577040e-02 
4.462882147989575e-02 
6.295732618092606e-02 
8.644035241970913e-02 
1.166164809306920e-01 
1.546690628394902e-01 
1.999554346680615e-01 
2.434683359132119e-01 
2.800846274146029e-01 
3.368595257878888e-01 
4.044418359833648e-01 
4.685002493634456e-01 
5.185062817085154e-01 
5.811314144990846e-01 
6.545700991450585e-01 
7.276588861478224e-01 
7.960626077582168e-01 
8.572037183403355e-01 
9.091330485015775e-01 
9.503131649503738e-01 
9.795718963793163e-01 
9.961006479199827e-01 



WEIGHTS 



9.878088201321919e-03 
1.109316819462674e-02 
1.313311581321880e-02 
1.624262442061470e-02 
2.065168462990214e-02 
2.657795406825320e-02 
3.399052299072427e-02 
4.208214612865170e-02 
4.732516974042797e-02 
3.618419415803922e-02 
4.547346840583578e-02 
6.463153575242817e-02 
6.859104457897808e-02 
5.589917935916451e-02 
5.199232318335285e-02 
7.089840644422261e-02 
7.427400331494240e-02 
7.125308736931726e-02 
6.513697474660338e-02 
5.682298546820264e-02 
4.678000924507099e-02 
3.538488886617123e-02 
2.299723483013955e-02 
9.993597414733579e-03 



24 point quadrature rule for integrals 
of the form J" 1 fix) + g(x) \og(x + x)dx, 
where 10~ 4 < x < 10" 3 



NODES 



2.625961371586153e-03 
6.309383772392260e-03 
1.073246133489697e-02 
1.645170499644402e-02 
2.433800511777796e-02 
3.582530925992294e-02 
5.315827372101662e-02 
7.917327903614484e-02 
1.162053707416708e-01 
1.648139164451449e-01 
2.231934088488800e-01 
2.864519293820641e-01 
3.466729491189400e-01 
4.076175535528108e-01 
4.800964107543535e-01 
5.594105009204460e-01 
6.395390292352857e-01 
7.167410782176877e-01 
7.882807127957939e-01 
8.519356675821297e-01 
9.058606177202579e-01 
9.485539755760567e-01 
9.788566874094059e-01 
9.959649506960162e-01 



WEIGHTS 



3.441901737135120e-03 
3.978799794732070e-03 
4.958449505644980e-03 
6.620822501994994e-03 
9.385496468197222e-03 
1.396512052439178e-02 
2.119383832447796e-02 
3.124989308824302e-02 
4.291481168916344e-02 
5.400832278279924e-02 
6.197424674301215e-02 
6.297221626131570e-02 
5.794981636764223e-02 
6.650501614478806e-02 
7.716379373230733e-02 
8.047814122759604e-02 
7.917822434973971e-02 
7.477646096014055e-02 
6.793424765652059e-02 
5.906852968947303e-02 
4.853108558910315e-02 
3.666228059710319e-02 
2.380850649522536e-02 
1.034186239262945e-02 



24 point quadrature rule for integrals 
of the form J" 1 fix) + g(x) log(x + x)dx, 
where 10 -5 < x < 10" 4 



NODES 


WEIGHTS 


7 


759451679242260e- 


04 


1 


049591733965263e-03 


1 


952854410117286e- 


03 


1 


314968855711329e-03 


3 


429053832116395e- 


03 


1 


651475072547296e-03 


5 


301128540262913e- 


03 


2 


135645684467029e-03 


7 


878118775220067e- 


03 


3 


165043382856636e-03 


1 


205537050949829e- 


02 


5 


479528688655274e-03 


1 


965871512055557e- 


02 


1 


028817002915096e-02 


3 


403328641997047e- 


02 


1 


923291785614007e-02 


5 


947430305925957e- 


02 


3 


212643438782854e-02 


9 


873500543531440e- 


02 


4 


638626850049229e-02 


1 


518862681939413e- 


01 


5 


960676923068444e-02 


2 


171724325134259e- 


01 


7 


052360405410943e-02 


2 


919941878735093e- 


01 


7 


863451090237836e-02 


3 


734637353255530e- 


01 


8 


381771698595157e-02 


4 


586710018443288e- 


01 


8 


612755554083525e-02 


5 


448057416999684e- 


01 


8 


569938467103264e-02 


6 


292158981939618e- 


01 


8 


271051499695768e-02 


7 


094415843889587e- 


01 


7 


736692567834522e-02 


7 


832417328632321e- 


01 


6 


990012937760461e-02 


8 


486194141302759e- 


01 


6 


056687669667680e-02 


9 


038469 149367938e- 


01 


4 


964868706783169e-02 


9 


474898150194623e- 


01 


3 


745026957972177e-02 


9 


784290662963747e- 


01 


2 


429741981889855e-02 


9 


958843370550371e- 


01 


1 


054906616108520e-02 
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24 point quadrature rule for integrals 






of the form J Q f(x) + g(x) \og(x + x)dx, 
where 1(T 6 < x < 1CT 5 








NODES 


WEIGHTS 


3 


126377187332637e-04 


4 


136479682893960e- 


04 


7 


671264269072 188e-04 


5 


068714387414649e- 


04 


1 


359575160544077e-03 


7 


008932527842778c- 


04 


2 


238313285727558e-03 


1 


110264922990352e- 


03 


3 


770276623583326e-03 


2 


120108385941761c- 


03 


7 


146583956092048e-03 


5 


249076343206215e- 


03 


1 


635515250548719e-02 


1 


450809938905405e- 


02 


3 


828062855101241e-02 


2 


987724029376343e- 


02 


7 


628984500206759e-02 


4 


593298717863718e- 


02 


1 


294255336121595e-01 


5 


98763447553802 le- 


02 


1 


949876755761554e-01 


7 


065953519392547e- 


02 


2 


693852297828856e-01 


7 


729918562776261e- 


02 


3 


469762441631538e-01 


7 


556635340171830e- 


02 


4 


122748928895491e-01 


5 


234123638339037e- 


02 


4 


662499202239145e-01 


6 


532130125393047e- 


02 


5 


421402737123784e-01 


8 


188272080198840e- 


02 


6 


248832413655412e-01 


8 


237354882288161e- 


02 


7 


053258496784840e-01 


7 


795795664563893e- 


02 


7 


798841313231049e-01 


7 


076514272025076e- 


02 


8 


461534275163378e-01 


6 


145788741452406e- 


02 


9 


022312524979976e-01 


5 


044339641339403e- 


02 


9 


465899812310277e-01 


3 


8078171 18430632e- 


02 


9 


780549563823810e-01 


2 


47154901 1101626e- 


02 


9 


958125149101927e-01 


1 


073289672726758e- 


02 



24 point quadrature rule for integrals 
of the form J Q f(x) + g(x) \og(x + x)dx, 
where 10~ 7 < x < 10" 6 



NODES 



1.019234906342863e-04 
2.506087227631447e-04 
4.461429005344285e-04 
7.422845421202523e-04 
1.289196091156456e-03 
2.739287668024851e-03 
9.075168969969708e-03 
2.968005234555358e-02 
6.781742979962609e-02 
1.217792474402805e-01 
1.886625378438471e-01 
2.650602155844836e-01 
3.465113608339080e-01 
4.178374197420536e-01 
4.597624982511183e-01 
5.348065111487157e-01 
6.194640153146728e-01 
7.013481004172354e-01 
7.770386175609082e-01 
8.442211768916794e-01 
9.010272836291835e-01 
9.459409782755001e-01 
9.777905486554876e-01 
9.957622871041650e-01 



WEIGHTS 



1.349775051746596e-04 
1.663411550150506e-04 
2.328782111562424e-04 
3.804721779784063e-04 
7.930350452911450e-04 
2.600694722423854e-03 
1.212249113599252e-02 
2.946708975720586e-02 
4.647771960691390e-02 
6.095376889009233e-02 
7.224844725827559e-02 
7.986429603884565e-02 
8.143206462900546e-02 
5.040529357007135e-02 
5.592137651001418e-02 
8.398073572656715e-02 
8.402586870225486e-02 
7.922223490159952e-02 
7.177919251691964e-02 
6.227551999401272e-02 
5.108407212719758e-02 
3.854783279333592e-02 
2.501496650831813e-02 
1.086176801402067e-02 





24 point quadrature rule for integrals 






of the form J* f{x) 


-f g(x) \og(x + x)dx, 






where 10 -8 


< 


x < 10~ 7 




NODES 


WEIGHTS 


3 


421721832247593e-05 


4 


559730842497453c- 


05 


8 


533906255442380e-05 


5 


840391255974745e- 


05 


1 


563524616155011e-04 


8 


761580900682040e- 


05 


2 


746612401575526e-04 


1 


617264666294872e- 


04 


5 


408643931265062e-04 


4 


433543035169213e- 


04 


1 


782382096488333e-03 


3 


116175111368442e- 


03 


1 


101243912052365e-02 


1 


655494413772595e- 


02 


3 


553172024884285e-02 


3 


242539256461602e- 


02 


7 


554170435463801e-02 


4 


734426463929677e- 


02 


1 


295711894941649e-01 


6 


032614603579952e- 


02 


1 


953213037793089e-01 


7 


069975187373848e- 


02 


2 


699680545714222e-01 


7 


806973621204365e- 


02 


3 


503697281371090e-01 


8 


216350598137868e- 


02 


4 


330838596494367e-01 


8 


261286657092808e- 


02 


5 


141801680435878e-01 


7 


883476216668445e- 


02 


5 


895097016206093e-01 


7 


157205125318401e- 


02 


6 


582708672338614e-01 


6 


703064468754417e- 


02 


7 


252543617887320e-01 


6 


706137273719630e- 


02 


7 


914154485613720e-01 


6 


4499841 16349734e- 


02 


8 


528383935857844e-01 


5 


775434959088197e- 


02 


9 


059696536862878e-01 


4 


812600239023880e- 


02 


9 


484664124578303e-01 


3 


661415869304224e- 


02 


9 


787863313133854e-01 


2 


386304203446463e- 


02 


9 


959482975155097e-01 


1 


038268695581411c- 


02 



24 point quadrature rule for integrals 
of the form J" 1 fix) + g(x) log(x + x)dx, 
where 10 -9 < x < 10" 8 



NODES 



6.538987938840374e-06 
2.613485075847413e-05 
5.664183720634991e-05 
1.179374114362569e-04 
3.299119431334128e-04 
3.626828607577001e-03 
2.265102906572155e-02 
5.896796231680340e-02 
1.092496277855923e-01 
1.666701689499393e-01 
2.196889385898800e-01 
2.770352260035617e-01 
3.483163928268329e-01 
4.153287664837260e-01 
4.695624219668608e-01 
5.421129318998841e-01 
6.238832212055707c-01 
7.041842972237081e-01 
7.788817007552110e-01 
8.453877637047045e-01 
9.017178251963006e-01 
9.462999385952402e-01 
9.779333485180249e-01 
9.957890687155009e-01 



WEIGHTS 



1.500332421093607e-05 
2.367234654253158e-05 
4.007286246706405e-05 
9.497743501485505e-05 
4.619067037944727e-04 
9.985382463808036e-03 
2.805741744607257e-02 
4.404106103008398e-02 
5.548413172821072e-02 
5.693235996372726e-02 
5.087307376046002e-02 
6.593729718379782e-02 
7.335680008972614e-02 
5.675029500743735e-02 
6.117926027541254e-02 
8.004805067067550e-02 
8.196991767042605e-02 
7.800219127200407e-02 
7.097175077519494e-02 
6.171193295041172e-02 
5.068671319716005e-02 
3.827738423897266e-02 
2.485063762733620e-02 
1. 0792849733295 16e-02 
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24 point quadrature rule for integrals 






24 point quadrature rule for integrals 


of the form f. f(x) 


+ g(x) \og(x + x)dx, 






of the form f. f(x) + g(x) log(x + x)dx, 


where 1CT 10 < x < 






where 10 -11 


< x < W- 1 " 


NODES 


WEIGHTS 


NODES 


WEIGHTS 


6.725520559705825e-06 


8 


128391913974039e-05 


2 


828736694877886e-08 


1 


665602686704325e-05 


6.986424152770461e-06 


-7.773900735768282e-05 


2 


302233157554212e-06 


2 


577419924039251e-06 


1.217363416714366e-05 


1 


287386499666193c- 


05 


5 


853587143444178e-06 


4 


957941112780975e-06 


2.677746219601529e-05 


1 


895577251914526e- 


05 


1 


451588770083244e-05 


1 


537074702915107e-05 


5.597036348896741e-05 


4 


732580352158076e- 


05 


9 


711965099273031e-05 


4 


640075239797995e-04 


2.729343280943077e-04 


9 


857909615386162e- 


04 


9 


004761967373848e-03 


1 


705687938176189e-02 


9.445526806263141e-03 


1 


756872897270054e- 


02 


3 


442077924035546e-02 


3 


349724914160473e-02 


3.556725025161542e-02 


3 


439422017906772e- 


02 


7 


543926781582543e-02 


4 


820210872119093e-02 


7.765556668177810e-02 


4 


944188361792970e- 


02 


1 


300373356318913e-01 


6 


054547286337976e-02 


1.336848150648662e-01 


6 


219733934997792e- 


02 


1 


955182772803384e-01 


6 


984354388121057e-02 


2.011576917683550e-01 


7 


228007436918939e- 


02 


2 


683608546664295e-01 


7 


498721497014774e-02 


2. 7727368543 14979e-01 


7 


944986391225688e- 


02 


3 


430029178740901e-01 


7 


240620145057083e-02 


3.590124362607926e-01 


8 


347646288178011e- 


02 


4 


085056107803621e-01 


5 


774925310174693e-02 


4.430074035214462e-01 


8 


380433020121207e- 


02 


4 


660198270439085e-01 


6 


238505554837956e-02 


5.247388219574510e-01 


7 


832768209682506e- 


02 


5 


336124745634699e-01 


6 


940394677081842e-02 


5.961053238782420e-01 


6 


300796225242940e- 


02 


5 


985245800106473e-01 


5 


910843483407385e-02 


6.547331131213409e-01 


5 


923406014585053e- 


02 


6 


564089719608276e-01 


6 


059752321454190e-02 


7.192258519628951e-01 


6 


834293563803810e- 


02 


7 


216666024232565e-01 


6 


823362237770209e-02 


7.874251789073102e-01 


6 


660337204499726e- 


02 


7 


893712241343741e-01 


6 


593839664071163e-02 


8.505852012775045e-01 


5 


911988751082552e- 


02 


8 


518883782001418e-01 


5 


853014420243146e-02 


9.047824617894323e-01 


4 


893575310568894e- 


02 


9 


055688088881344e-01 


4 


849217100974983e-02 


9.479045131744448e-01 


3 


708256438629509e- 


02 


9 


483163097840529e-01 


3 


677417821170115e-02 


9.785770588866582e-01 


2 


411463784693618e- 


02 


9 


787413692715607e-01 


2 


392585642844202e-02 


9.959104692340199e-01 


1 


048087156697020c- 


02 


9 


959413203611228e-01 


1 


040149939671874e-02 



24 point quadrature rule for integrals 






24 point quadrature rule for integrals 


of the form J Q f(x) 


+ < 


j(x) \og(x + x)dx, 






of the form J f(x) 


+ g(x) log(x + x)dx, 


where 10" 12 


<x< 10~ n 






where 10" 13 


< 


x < 10" 12 


NODES 


WEIGHTS 


NODES 


WEIGHTS 


6.147063879573664e-07 


8 


763741095000331e-07 


I 


523740015216508e-08 


4 


418138082366788e-07 


2.102921984985835e-06 


1.784696796288373e-05 


4 


281855233588279e-07 


4 


389108058643120e-07 


2.188366117432289e-06 


-1.795398395983826e-05 


1 


036900153156159e-06 


9 


539585150737866e-07 


3.482602942694880e-06 


5 


117514567175025e- 


06 


7 


825849325746907e-06 


5 


823980947200484e-05 


2.768001888608636e-05 


1 


698863549284390e- 


04 


8 


617419723953112e-03 


1 


634464263521301e-02 


8.942779215792784e-03 


1 


701975216672032e- 


02 


3 


268881 163637599e-02 


3 


129682188728318e-02 


3.432218364237253e-02 


3 


346025972593909e- 


02 


6 


988441391437043e-02 


4 


212468617589480e-02 


7.530931328026620e-02 


4 


817949622196712e- 


02 


1 


142202307676442e-01 


4 


505120897719191e-02 


1.298983048592572e-01 


6 


055152664710045e- 


02 


1 


596471081833281e-01 


4 


769069780026684e-02 


1.954020797117703e-01 


6 


988313730886592e- 


02 


2 


135336418959620e-01 


6 


038503382768951e-02 


2.682970870436427e-01 


7 


504602275463067e- 


02 


2 


781100275296151e-01 


6 


695343672694180e-02 


3.429540704041702e-01 


7 


230942674874111e- 


02 


3 


433392803364457e-01 


6 


163298712826237e-02 


4.080399755202422e-01 


5 


705952259766429e- 


02 


4 


019960595528027e-01 


5 


877742624357513e-02 


4.652562798154792e-01 


6 


265021180818162e- 


02 


4 


656415679416787e-01 


6 


800053637773440e-02 


5.333220999210325e-01 


6 


993669694523695e- 


02 


5 


334880548894250e-01 


6 


516918103589647e-02 


5.986982369433125e-01 


5 


937130986945129e- 


02 


5 


943298528903542e-01 


5 


853785375926075e-02 


6.564773600603511e-01 


6 


026572020863567e- 


02 


6 


562968737815924e-01 


6 


639396325654251e-02 


7.215159032030418e-01 


6 


815292696374753e- 


02 


7 


25034334460 1498e-01 


6 


948738324081696e-02 


7.892098210760941e-01 


6 


596804590657802e- 


02 


7 


928820737781136e-01 


6 


538801703374268e-02 


8.517672777806986e-01 


5 


857483758149194e- 


02 


8 


546103048745466e-01 


5 


761503751629250e-02 


9.054906995605498e-01 


4 


853209199396977e- 


02 


9 


073762310762705e-01 


4 


761344859555310e-02 


9.482736017320823e-01 


3 


680469214176019e- 


02 


9 


493253659835347e-01 


3 


607033097268266e-02 


9.787238593479314e-01 


2 


394561701705853e- 


02 


9 


791606801267259e-01 


2 


345690720840071e-02 


9.959379852805677e-01 


1 


041005152890511e- 


02 


9 


960217573957566e-01 


1 


019557402722854e-02 



24 point quadrature rule for integrals 
of the form J Q f(x) + g(x) \og(x + x)dx, 
where 1(T 14 < x < 1(T 13 



NODES 

6.025980282801020e-08 
6.411245262925473e-08 
1.862815529429129e-07 
2.029190208906422e-06 
8.902881307076499e-03 
3.420089035164912e-02 
7.508687525931594c-02 
1.295858123029775e-01 
1.950409815188335e-01 
2.679751967812604e-01 
3.428525062164689e-01 
4.080941369413548e-01 
4.64664451 1900009e-01 
5.328071517215501e-01 
5.978508749698001e-01 
6.521214523350964e-01 
7.134921670665336e-01 
7.679317896479284e-01 
8.029718487208403e-01 
8.551101435866935e-01 
9.067319102017767e-01 
9.487765213293372e-01 
9.788979796532736e-01 
9.959684838634199e-01 



WEIGHTS 
9.079353616441234e-07 
-8.390389042773805e-07 
2.782460677485016e-07 
1.821115881362725e-05 
1.695809650660321e-02 
3.336370146025145e-02 
4.807898681796971e-02 
6.047672723211479e-02 
6.986774906175534e-02 
7.515608233194288e-02 
7.264249904037610e-02 
5.672507168477261e-02 
6.220316364524964e-02 
7.032362652293805e-02 
5.742730804758014e-02 
5.644075454541 152e-02 
6.318643666150391e-02 
3.945995610428228e-02 
4.324200884758527e-02 
5.478223695609097e-02 
4.740856250832772e-02 
3.633314063504751e-02 
2.372788917088821e-02 
1.033036588606145e-02 



